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Abstract 


We consider a class of singularly perturbed semilinear three-point boundary 
value problems. An accelerated uniformly convergent numerical method is 
constructed via the exponential fitted operator method using Richardson 
extrapolation techniques to solve the problem. To treat the semilinear 
term, we use quasi-linearization techniques. The numerical results are tab- 
ulated in terms of maximum absolute errors and rate of convergence, and 
it is observed that the present method is more accurate and ¢-uniformly 
convergent for h > ¢, where the classical numerical methods fail to give a 
good result. It also improves the results of the methods existing in the lit- 
erature. The method is shown to be second-order convergent independent 
of perturbation parameter e. 
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1 Introduction 


Singularly perturbed differential equations are typically characterized by the 
presence of a small positive parameter ¢ multiplying some or all of the highest 
order terms in differential equations. Such types of problems arise frequently 
in mathematical models of different areas of physics, chemistry, biology, en- 
gineering science, economics, and even sociology. 


Boundary value problems including nonlocal conditions, which connect 
the values of the unknown solution at the boundary with values in the in- 
terior, are known as nonlocal boundary value problems. The study of this 
kind of problem was initiated by Il’in and Miseev [18, 19], motivated by the 
work of Bitsadze and Samarskii [4] on nonlocal linear elliptic boundary value 
problems. These problems have been used to represent mathematical mod- 
els of a large number of phenomena, such as problems of semiconductors 
in electronics, the vibrations of a guy wire of a uniform cross-section, heat 
transfer problems, problems of hydromechanics, catalytic processes in chem- 
istry and biology, the diffusion-drift model of semiconducting devices, and 
some other physical phenomena; see [1, 17, 25]. The existence and unique- 
ness of the solutions to nonlocal boundary value problems have been studied 
by many authors [3, 20]. Some approaches for the numerical solution of sin- 
gularly perturbed nonlocal boundary value problems have been proposed in 
(2, 6, 7, 11, 12, 13, 16, 21, 26]. Uniformly convergent numerical methods 
of order second and high for solving different singularly perturbed problems 
have been studied in [5, 9, 10, 14, 23, 27]. What makes the problem under 
consideration challenging is that, it contains the perturbation parameter in 
both diffusion and convection terms, that it contains semilinear source term, 
and that it has three-point nonlocal boundary condition. To the best of our 
knowledge, the problem under consideration has not been done using the ex- 
ponentially fitted operator method. Motivated by the paper [8], we develop a 
uniformly convergent numerical method for solving the singularly perturbed 
problem under consideration. 


2 Definition of the problem 


Consider the following singularly perturbed problem of the form 


Ly := €?y"(x) + ea(x)y'(x) — f(x, y(a)) =0, O<2<l, (1) 

with the given conditions 
y(0) = A, (2) 
Loy := y(l) — d(y()) = 0, O<h <I, (3) 
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where ¢,0 <¢ <1 is the perturbation parameter, A is a given constant, and 
the functions a(x) > a > 0 and f(x,y) are sufficiently smooth on (0, /] and 
[0,2] x R, respectively. Moreover 0 < b < ae <6 <oo and || <k<1. 

The solution of y(x) of (1)-(8) has boundary layers at x = 0 and x = 1 
for € near 0. 

Equations of this type arise in mathematical problems in many areas of 
mechanics and physics. Among these are the Navier-Stokes equations of 
fluid flow at high Reynolds number, mathematical models of liquid crystal 
materials and chemical reactions, shear in second-order fluids, control theory, 
electrical networks, and other physical models [23, 24]. 

To obtain the numerical solution to (1)-(3), the well-known Newton’s 
quasi-linearization techniques were used. This technique allows us to linearize 
the semilinear problem into a linear problem, whose solution y“) (x) with a 
proper initial guess y) will converge to the original solution y(2). 

So, after applying quasi-linearization technique for f(x, y(a)), we rewrite 
(1)-(3) in the form 


Ly = ey" («) + ea(a)y'(x) —W(a)y(a) = R(x), rE N=(0,0, (4) 


with the given conditions 
y(0) = A, (5) 


Loy := y() — o(y(h)) =0, 0< <1, (6) 


where b(x) = $4(a,ny(a)),0<m <1, and R(x) = f(x, ny(x)) — ny(w)b(2). 

In this article, we analyze an exponentially fitted operator scheme with 
Richardson extrapolation techniques on uniform mesh for the numerical so- 
lution of (4)-(6). Uniform convergence is proved in the discrete maximum 
norm. Finally, we formulate the algorithm for solving the discrete problem 
and give the illustrative numerical results. 


3 Properties of continuous solution 


The following lemmas are necessary for the existence and uniqueness of the 
solution and for the problem to be well-posed. 


Lemma 1. Let y(x) be the solution of (1)-(3). If a € C'[0,l], then the 
following inequality 

IIylloo S Co (7) 
holds, where Cp = (1—)71(|A]+|B|+07"||R||0), B = $(0,) R(x) = f(a, 0), 
and ||y|lo0 = maxo,y |y(*)]- 


Proof. For proof refer [8]. 


Lemma 2. Let y- be the solution of (P-). Then, for k = 0,1, 2,3, 4, 
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ja) Is © (1 +e°* (exp (=2") ex »(= “—")))), x € (0,0, 
of 
Oy 


(8) 
holds fore the solution y(a) provided that —ea'(x) > b and gf <C 
for x € [0,l] and |y| < C, where co = $[,/a?(0) + 4b + a(0)] and c = 
5[\/a?(1) + 4b — a(l)}. 


Proof. For proof refer [8]. 


4 Formulation of numerical scheme 


We subdivide the domain 2 = [0, 1] into N equal number of subintervals, each 
of length h. In this article, we develop the fitted operator finite difference 
method to find a numerical solution to the problem (4). We use the theory 
in the asymptotic method for developing the exponential factor. 

In order to evaluate the fitting factor, we divide both sides of (4) by «, 
and we obtain 


Ly(x) = ey" (x) + a(x)y' + p(x)y(z) = g(x), O<a<Z, 
= 44, (9) 
(0) =90, 0<4 <8, 


where p(«) =~", g(x) = 2, and 6 = g(y(é1)). 


4.1 Left end boundary layer problem 


The boundary layer occurs on the left side of the domain, that is, near x = 0. 
To find the numerical solution of (9), we use the theory applied in the asymp- 
totic method for solving singularly perturbed boundary value problems. From 
the theory of singular perturbations given by O’Malley [24] and using Tay- 
lor’s series expansion for a(x) about x = 0 and restriction to their first terms, 
we get the asymptotic solution as 


y(2) = yo(e) + D(A — yo(o)ye~F AP 8B) +O), (10) 


a(a 


where Yyo(x) is the solution of the reduced problem (obtained by setting ¢ = 0) 
of (9) given by a(x)yo(x) + p(x)yo(2) = g(a), with y(0) = A and yo(l) = 8, 
where 0 = y(y(I1)) 

By taking Taylor series expansion for a(x) and p(a) about the point “0” 
and simplifying them, we obtain 
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a(0) _ p(O) 


y(2) = yo(@) + (A yo(0))e7 AP 83)* + Of). (11) 


Now we divide the interval [0, 1] into N equal subinterval of mesh size h = + 


so that 2; = th, i=0,1,2...,N. From (11), we have 


a(0) _ p(0)) « 
ylih) = yo(ih) + (A —yo(0))e7 SH) 4. O(e), 
Taking limit as h — 0 on both sides, we have 
. : —(2@._ 20) ip 
lim y(ih) = yo(0) + (A= yo(o)e" E84 Of). (12) 


Also, we discretized 1341 = (+ 1)h =ih+h from (11), we obtain 


a(0) _ p(0) 


ylih + h) = yo(ih + h) + (A— yo(0))e7 PAH)" + O(6), 


Taking limit as h + 0 on both sides, we have 


a(0)_ p 


lim y(ih +h) = yo(0) + (4 — yo(0))e“ EP 


ite + O(6), (13) 
Similarly for x;_; = (i — 1)h = th — h from (11), we have 


a(0) _ 


ylih —h) = yo(éh — h) + (A — yo(O))eW FPF)" 4 OL), 


Taking limit as h - 0 on both sides, we have 


lim y(ih — h) = yo(0) + (A—yo(oyJe CP BB) °-# 4 Of). (14) 


To handle the effect of the perturbation parameter, artificial viscosity (expo- 
nentially fitting factor o(p)) is multiplied on the term containing the pertur- 
bation parameter as 


ea(p)y; + a(xi)y; + p(ai)ys = g(2), (15) 


with boundary conditions y(0) = A and y(@)=6. Next, we consider the 


=N 
difference approximation of (9) on a uniform grid Q = {a;} and denote 
h= Li4+1 — Vj. 


Using the Taylor series expansion for y(a) about the point x;, this can be 
written as 


h2 

yah 4 =e hy; Tr ai a O(h*), (16) 
h? 3 

yep th) > yin =yi thy, t ore + O(h”). (17) 
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For any mesh function y;, from (16) and (17), we obtain the following finite 
difference operator: 


D~y; = an ee 
+,,. — Yitl Vi 
D Yi 3 h ’ (18) 
D°y; = Yip i 


2 
+p-,,. — yitir=2yityi-1 


? 


Then, applying the central and forward finite difference formula for the second 
and the first derivative, respectively, on (9), we get 


eo(p)(D* D~y(a;)) + a(ai)(DTy(ai)) + p(wi)y(@i) = g(ai). (19) 
Using the difference operator in (19), we have 
Le ys = 9, (20) 


with boundary conditions y(0) = A and y(£) = 0. 


From (19), we have 


i-1 — 2yi + yi +1 — Yi 
co(p)(4 . Ht) a(n) (HA 8) a pledyn = ole) 1 


IA 
IA 


wl eS 


where p = A, 


Multiplying (21) by h, considering h small, and truncating the term 
(9(xi) — p(wi)yi)h, we get 


7 (uit ~ yi + yisr) + (xi) (yis1 — yi) =A(g(ai) — pigs). (22) 


By evaluating the limit of (22) as h approaches to zero and limp+o h(gi — 
piyi) = 0, we get 


limp+o(Yit1 — Yi) 
: : 23 
( limp so (Yin1 — 2y; + yi-1) 8) 


C= 


Substituting (12), (13), and (14) into (23) for solving the fitting parameter 
o and simplifying it, we obtain 


2 
a? (0)—eb(0) 
m0 P) 


4 | sink? (“Soy I 
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4.2 Right end boundary layer problem 


The boundary layer occurs on the right side of the domain, that is, near x = 1. 
Assume that @” denotes the partition of [0,4] into N subintervals such that 
0=20, %1,...,2n =1 with a =ih, h= % i=0,1,2,...,N. To find the 
numerical solution of (9), we use the theory applied in the asymptotic method 
for solving singularly perturbed boundary value problems. Assume that the 
problem has the right layer. From the theory of singular perturbations given 
by O’Malley [24] and using Taylor’s series expansion for a(x) about 2 = 1 
and restriction to their first terms, we get the asymptotic solution as 


a(1) a i ( 
= —_——_ A) _— 1 xc € 
ula) = wo(2) + 2 (0 wo(1))¢ 
where yo(z) is the solution of the reduced problem (obtained by setting « = 0) 
of (9) given by a(x) y6(x)+p(x)yo(x) = g(a), with yo(1) = 6. By taking Taylor 
series expansion for a(x) and b(x) about the point “1” and simplifying them, 
we obtain 


a(x) _ p(x) 


3) de Ole), (25) 


a(1)(1—2) 


u(x) = yola) + 0 yo(t)en" CPE) + Of). (26) 


Now we divide the interval [0, 1] into N equal subintervals of mesh size h = x 
so that x; = th, i=0,1,2...,N. From (26), we have 


lim y(ih) = yo(0) + (6 yo(t)ye—"M2) 406). 27) 
Also, we discretization x41 = (¢+1)h = ih +h From (26), we have 
lim y(ih + A) = yo(0) + 8 yolte@ME-) + Of). (28) 
Similarly for 7;-1 = (¢- 1)h = ih —h, we have 


lim y(ih — h) = yo(0) + (8 wo(tye2ME—) 4 (6). (29) 


Now we consider the second-order upwind finite difference scheme from (15). 
We have 


Yi-1 — Yi + Yi Yi — Yi- 
ca(p)( : he +) a a (4) + PiYi = Gi- (30) 


Multiplying (30) by h, considering h small, and truncating the term (g(z;) — 
p(xi)yi)h, we get 


7 int — 2y; + yira) +a:(ys — yi-1) = h(gi — Diy). (31) 
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By evaluating the limit of (31) as h approaches to zero and using limp_4o h(gi— 
piyi) = 0, we get 


limp—+o0(yi — Yi-1) 
limp+o(Yi—1 — 2ys + Yori) 


(32) 


g=-—a(l 
Substituting (27), (28), and (29) into (32) for solving the fitting parameter 
o and rearranging it, we obtain 
ear 5 _ 1) 


o> = pal) | (e- (33) 


4 | sink? (“Saya 


Case (1): For the left layer, from the discrete form of (4) on the domain 
[1, Ny, we have 
ery! + EQiy; = biyG = R;. (34) 


Using the central and forward finite difference formula on (34), we have 


i-1 — 2Yi + Yi +1 — Yi 
2o(e)(4 — HA) 4 oq (He 8) — bays = Be 


h2 


Hence, the required finite difference scheme becomes 


<*01(p) —2e*o1(p) a eo1(p) | ai 
( h2 Jaa + h2 E is b; Yi } he FE h Yi41 = R;. 


(35) 
The numerical scheme in (35) can be written as a three-term recurrence 
relation by 


Byithyt+ Gyn =i, 1sis< > (36) 
where 
€*01(p) —2e701(p) Qj 
EB = F; = bj ) 
i= (Gee), n= (Gre 
€*01(p) a; 
= gar Mls A; = R; 
a= (“rb +eF) 


Case (2): For the right layer, from the discrete form of (4) on the domain 
(, N — 1], we have 
e*yy + cary), — biys = Rj. (37) 


Using the central and backward finite difference formula on (37), we have 


i-1 — 2yi + yi i — Yi- 
oa(o)( 4 ; a) reai( 4 7 :) biys = Ri. 
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Hence, the required finite difference scheme becomes 


€*02(p) __a(zi) —2e*02(p) a €*02(p) 
OP encase 
(38) 


The numerical scheme in (38) can be written as a three-term recurrence 
relation by 


N 
Eiyyi-1 + Fiys + Giyinn = Ai, a i<N-I, (39) 


where 


<*02(p) 


Next, in order to treat the boundary conditions, the following equations are 
obtained: For i = 1, (36) becomes 


Fyy. + Giy2 => Ay = BE, A. (40) 
For i = N — 1, (36) becomes 


En_-iyn—-2 + Fn-i1yn-1+ Gn-iyn = Hn-1. (41) 


From (6) and (41), we obtain 


En_-iyn-2 + Fn—-iyn-1+ Gn 1O(yx) = Hy-1. (42) 


Therefore, the problem in (1) with given boundary conditions in (2)—(3), can 
be solved using the scheme (36), (39), (40), and (42), which forms an N x N 
system of algebraic equations. 


5 Uniform convergence analysis 
In this section, we need to show the discrete scheme in (36) and (39) satisfies 
the discrete minimum principle and uniform convergence. 


Lemma 3 (Discrete Minimum Principle). Let Y; be any mesh function that 
satisfies Yy > 0,Yy > 0 and L.Y; < 0, i =1,2,3,...,N—1. Then Y; >0, 
for i=0,1,2,...,N. 


Proof. The proof is by contradiction. Let 7 be such that Y; = min, Y;, and 
suppose that Y; < 0. Clearly, 7 ¢ {0,N}. Yj41 — Yj > 0 and Y; — Yj-1 < 0. 
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Therefore, 
iy ae Vou 
LY; -<( a | +4; (A, 2) + pi¥), 

a i i Qj i! 

=F (eat —2Y74- ei) 4 5 (Yj41 — Yj) + Dj, (43) 
E a; 

= 72 ((V41 — ¥5) — (¥ — ¥j-1)) + 7, Wi — ¥;) + p5¥3, 

>0 


* 


where the strict inequality holds if Y;,,; — Y; > 0. This is a contradiction, 
and therefore Y; > 0. Since j is arbitrary, we have Y; > 0, for 7 = 
0,1,2,...,N. 


We proved above that the discrete operator Le satisfies the minimum 
principle. Next, we analyze the uniform convergence analysis. 


Theorem 1. Let y(x;) and Y; be, respectively, the exact solution of (9) 
and numerical solutions of (20). Then for sufficiently large N, the following 
parameter uniform error estimate holds: 


IL’ (y(a;,) — Y;)| < Ch (1 +e max (exp (=S*)))) . (44) 


<N-1 
Proof. Let us consider the local truncation error defined as 


L™ (y(ai) — Yi) = e02(p)(y" (ai) — DE D~y(ai)) + a(ai)(y'(@i) — 2 
a? (a) |b) ») : ( ) 


@) =e 
— , since p= +. 
sinh? (G7?) p)) | E 


~( 
where €09(p) = 24 c 


For fixed N = i, taking the limit for ¢ > 0, we obtain 


a?(1)—eb(1) 


—(—aay ~ P) 
. . pa(l e a(1) —1 
lim eo(p) = lim iy c 3, a2 (nob) ) = Ch, 
e>0 e30 4 sinh? |( Zatt) ) 


where C is constant independent of h and «. 
Using Taylor series expansion, the bound for y(a;_-1) and y(aj+41) at x; is 
obtained as 


y(i—1) = y(ws) — hy (aa) + Sry" (@e) — Fry (wa) + ry (a) + OCF), 
y (wisi) = y(as) + hy! (wi) + Fy!" (wi) + Fy (wa) + Fy (ai) + O(R?). 


We obtain the bound for 


\D+ D-y(2)| < Cly" (aah 
{ ly" (as) - D+ D-y(,)| < Ch?|y (2,)}. ai 
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Similarly, for the first derivative term, we have 
ly'(xi) — DT y(2i)| < Chly"(ai)I, (47) 


where ly*) (x,)| = SUP, € (29,2) ly (zl, k = 2,3,4. 
Using the bounds in (46) and (47), we obtain 


IZ” (y(ai) — ¥i) < Ch? ly (wi)| + Cly”@a)- 


Now, using the bounds for the derivatives of the solution in Lemma (2), we 
have 


IL" (y(w;) — ¥i)| < Ch (1 +e" exp (S) 


& 


+Ch (1 +e exp (=) 


<Ch (1 +e-* max exp (=) 
1<i<N-1 E 


since e~4 > €7?, 


Lemma 4. For a fixed mesh and for ¢ > 0, it holds 


ary 
exp (= ) 
lim max —m——~*=0, m=1,2,3,..., 
e+01<i<N-1 em 


(=4==2) 

exp | — 

ina age pO, 
E70 1<i< N-1 em 


Proof. Consider the partition [0,1] := {0 = 2% < a1 <-++-<ay_-1<ay = 
1} for the interior grid points, we have 


1<i<N-1 em ~ Ene er 
—a(1— 2;) —a(1— y_-1) —a 
exp exp exp | —— 
€ E 
max < = 
1<i<N-1 em ai em 


as £1; =1—2y_,= h. The repeated application of L’Hospital’s rule gives 


exp | —— rs 
E r m! 


lim ——~—~+ = li —— = li —____________ = 9) 
e—30 em (2 exp(ahr) tog (ah)™ exp(ahr) 
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Theorem 2. Let y(x;) and Y; be the exact solution of (9) and numerical 
solutions of (20), respectively. Then, the following error bound holds 


sup |(y(ai) — Yi)| < Ch. (48) 
O0<e<<l 


Proof. By substituting the results in Lemma (4) into Theorem (1) and ap- 
plying the discrete minimum principle, we obtain the required bound. 


Richardson extrapolation 
Here, we apply the Richardson extrapolation technique to accelerate the rate 
of convergence of the proposed scheme. Richardson Extrapolation is a conver- 
gence acceleration technique, which involves a combination of two computed 
approximations of solution. From (48), we have 


y(zi) — Yi < Ch, (49) 


where y(z;) and Y; are exact and approximate solutions, respectively, and C’ 
is constant independent of ¢ and h. Let 2?" be the mesh obtained by bisecting 
each mesh interval in Q” and denote the approximation of the solution on 
0?" by Y2". From (49), we have 


y(xi) -Yi 2 Ch+R", 2 €ON, (50) 


So, this works for any B # 0 gives 
h 
y(2,) = yer ~ CS 4 a Li € o-, (51) 


where the remainders, R” and R?" are O(h?). Equations (50) and (51) lead 
to 
y(wi) — (2¥7" — ¥;) & Ch’, 


which gives that 
bial = oy" _ Y; (52) 


is also an approximate solution. The total truncation error for the approxi- 
mate solution in (52) becomes 


sup |(y(x:) — Yi)| < Ch’. (53) 
O<e<<l 
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6 Numerical example and results 


To validate the established theoretical results, we perform numerical experi- 
ments using the model problems of the form in (1)—(3). 


Example 1. Consider the model singularly perturbed boundary value prob- 
lem 


ey" (x) + e(1— Sy (2) —(1+2? + 2y—exp(-y))=0, O0<2<1, 


subject to the boundary conditions 


Having y; = yp (the approximated solution obtained via fitted operator 
finite difference method) for different values of h and ¢, the maximum errors. 
Since the exact solution is not available, the maximum errors (denoted by 
EP) are evaluated using the double mesh principle for fitted operator finite 
difference methods using formula 


N 


h 
e I. 


h 2 
= max = : 
pmax ly; Y25 


Furthermore, we will tabulate the ¢-uniform error 


EN = max a. 
O<e<l 


The numerical rate of convergence and the e-uniform rate of convergence are 
computed using the formulas 


ny _ log(E%) — log(E?%) 


ss log(2) 


E N= 16 N = 32 N = 64 N = 128 N = 256 
10-4 —-3.3238e-03 8.4445e-04 2.1268e-04 5.3355e-05 1.3361e-05 
10-8 —3.3238e-03  8.4445e-04 2.1268e-04 5.3355e-05 —1.3361e-05 
10~'2  3.3238e-03  8.44450-04 2.1268e-04 5.3355e-05 —1.3361e-05 
10~'© = 3.3238e-03  8.44450-04 2.1268e-04 5.3355e-05 —1.3361e-05 
10~7°  3.3238e-03  8.4445¢-04 2.1268e-04 5.3355e-05 —1.3361e-05 
EN — 3.3238e-03  8.4445e-04  2.1268e-04 5.3355e-05 —-1.3361e-05 
RY 1.9767 1.9893 1.9950 1.9976 


Table 1: Maximum absolute errors and rates of convergence for Example 1. 
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N=16 N = 32 N = 64 N = 128 N = 256 
Present 
EN 3.3238e-03 8.4445e-04 2.1268c-04 5.3355e-05  1.3361e-05 
RN 1.9767 1.9893 1.9950 1.9976 
Method in [8] 
EN 0.01388073  0.00955129 0.00550127 0.00314426 0.00176922 
RN 0.69 0.75 0.81 0.91 


Table 2: Comparison of maximum absolute errors and rate of convergence for Example 
1 at different number of mesh points. 


N= 16 N = 32 N = 64 N = 128 N = 256 


After 
EN 3.3238e-03  8.4445e-04 2.1268¢-04 5.3355e-05 1.3361e-05 
RN 1.9767 1.9893 1.9950 1.9976 

Before 
EN 6.1832e-02 3.2578e-02 1.6711e-02 8.4619e-03 4.2576e-03 
RN 0.9245 0.9631 0.9817 0.9909 


Table 3: Maximum absolute errors and rate of convergence before and after Richardson 
extrapolation for Example 1 at different number of mesh points. 


—— o(t/n) 


Maximum pointwise error 
3 

Maximum pointwise error 
3 


10! 10° 10° 10' 10° 10° 
N-Number of mesh points N- Number of mesh points 


Figure 1: e-uniform convergence of exponentially fitted operator method in log-log 
scale before and after Richardson Extrapolation, respectively. 


7 Discussion and conclusion 


This study introduced an accelerated exponentially fitted finite difference nu- 
merical method for solving singularly perturbed semilinear differential equa- 
tions with nonlocal boundary conditions. The behavior of the continuous 
solution of the problem was studied and shown that it satisfies the continu- 
ous stability estimate and the derivatives of the solution are also bounded. 
The numerical scheme was developed on the uniform mesh using exponential 
fitted operator method in the given differential equation. The stability of the 
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Figure 2: The e-uniform convergence of the method in log-log scale before and after 
extrapolation at the same plot for Example 1. 


developed numerical method was established, and its uniform convergence 
was proved. To validate the applicability of the method, a model problem 
was considered for numerical experimentation for different values of the per- 
turbation parameter and mesh points. The numerical results were tabulated 
in terms of maximum absolute errors, numerical rate of convergence, and 
uniform errors (see Table 1) and to show the performance of our scheme 
after Richardson extrapolation, we compared both the maximum absolute 
error and the rate of convergence of the scheme before and after Richardson 
extrapolation (see Table 3). To visualize our result more, we plotted the log- 
log plot of maximum absolute error before and after Richardson in different 
plane and in the same plane in Figures 1 and 2, respectively. The method 
was shown to be ¢-uniformly convergent with the order of convergence O(h?). 
The performance of the proposed scheme was investigated by comparing it 
with prior study (see Table 2). The proposed method gave more accurate, 
stable, and ¢-uniform numerical results. 
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